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Abstract 


We report a gradual brightening y-ray source, 4FGL J1718.5+4237, in 0.1-500.0 GeV, which may be associated 
with a blazar NVSS J171822+423948 with a redshift of ~2.7. We analyzed 15.25 yr of y-ray data recorded by the 
Large Area Telescope on board the Fermi Gamma-ray Space Telescope and detected significant y-ray emissions in 
the direction of the blazar with a test statistic (TS) of ~135. Based on timing analysis using a | yr time bin, we 
have observed a gradual brightening in y-ray emissions from the target. In our analysis, we categorize them into 
two states: Quiet (TS ~ 0) and Loud (TS ~ 226) states, with the distinction occurring in 2016 August (MJD 
57602.69). From the Quiet state to the brightest period (the last data point), the y-ray flux in 0.1-500.0 GeV 
increased by more than 12-fold from <0.2 x10 *photonscm™'s~' to 2.6 x 107° photons cm™' s™'. 
Additionally, we studied the spectral properties in detail for the Loud state and the overall data. While no 
significant variation was detected, both exhibited a spectral index I’ of ~2.6. The origin of the brightening y-ray 
emissions from the target is not yet clear. Future long-term multi-wavelength observations and studies will provide 
insight into the astrophysical mechanisms of the target. 
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1. Introduction 


As the most luminous sources across the entire sky, active 
galactic nuclei (AGNs) are energetic astrophysical sources 
powered by accretion onto supermassive black holes in 
galaxies, and present unique observational signatures in which 
multi-wavelength electromagnetic spectral energy distributions 
(SEDs) extend from MHz radio frequencies to TeV y-ray 
energies (Urry & Padovani 1995; Ulrich et al. 1997). The 
broadband SED from an AGN reveals a universal double broad 
hump structure in a logy—logvF,, diagram, with the low-energy 
hump covering the radio to ultraviolet (UV) and X-ray 
frequencies, mostly due to synchrotron emission from electrons 
in a relativistic jet, while the origin of the high-energy hump 
from X-rays to y-rays is highly debated. This component could 
be produced via inverse Compton (IC) scattering in the 
conventional leptonic emission scenario (e.g., Ghisellini et al. 
1985; Begelman & Sikora 1987; Maraschi et al. 1992; Dermer 
& Schlickeiser 1993; Sikora et al. 1994; Bloom & Marscher 
1996; Btazejowski et al. 2000; Finke et al. 2008; Ghisellini & 
Tavecchio 2009; Abdo et al. 2010a; Kang et al. 2014; Yan 
et al. 2014) or could be either via proton synchrotron emission 
inside the jet (e.g., Mücke & Protheroe 2001) or secondary 
particles in electromagnetic cascades initiated by pion decay in 
the hadronic emission scenario (e.g., Mannheim & Biermann 
1992; Aharonian 2000; Miicke & Protheroe 2001; Miicke et al. 


2003; Böttcher et al. 2013; Petropoulou & Mastichiadis 2015; 
Gasparyan et al. 2022). It is general believed that AGNs are 
classified as different subclasses based on their observational 
characteristics and physical properties. 

The Large Area Telescope on board the Fermi Gamma-ray 
Space Telescope (Fermi-LAT; Atwood et al. 2009) monitors 
the entire high-energy y-ray sky every three hours, which 
affords us an ideal opportunity to investigate the y-ray 
emissions from AGNs, unbiased by activity state or AGN 
sub-class. This ability has revealed the AGNs are dominant 
sources in the extragalactic y-ray sky. In the Fourth Fermi-LAT 
Source catalog (4FGL-DR4; Ballet et al. 2023), about 4015 
AGNs were published with energy >100 MeV, accounting for 
more than 55% of all Fermi-LAT sources. Among them, 
blazars, a subclass of AGNs with powerful relativistic jets 
oriented at small angles (<10°) with respect to the line of sight 
to the observer (Urry & Padovani 1995), are the most abundant 
proportion of known y-ray AGNs (Ballet et al. 2023). Blazars 
are typically subdivided into flat spectrum radio quasars 
(FSRQs), which display strong emission line components in 
the optical spectrum, and BL Lacertae objects (BL Lacs) that 
have weak or no emission lines (Urry & Padovani 1995). 
Among Fermi blazars, BL Lacs and FSRQs account for ~40% 
and ~20% of the sources, respectively, while blazar candidates 
of unknown types are about 35% of all Fermi blazars (3LAC; 
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Ajello et al. 2022). So far, a handful of y-ray blazars have been 
even found at very high-redshifts (z > 2) (e.g., Romani et al. 
2004; Abdo et al. 2010b; Ackermann et al. 2015, 2017; Paliya 
2015; Ajello et al. 2016, 2022; Paliya et al. 2016, 2019; Kaur 
et al. 2017, 2018; Liao et al. 2018; Sahakyan et al. 2020), and 
one of the major characteristics of these sources is that they 
have softer y-ray spectra. Here we present the y-ray emissions 
of 4FGL J1718.5+4237 from 4FGL-DR4 (Ballet et al. 2023), 
which may be associated with a blazar, NVSS J171822 
+423948, with a redshift of ~2.7. 

Zhang et al. (2023) carried out an analysis for y-rays from 
the globular cluster NGC 6341 with Fermi-LAT, and found a 
significant y-ray excess presented with a maximum likelihood 
test statistic (TS) of ~40 (see their paper for more details). 
Motivated by their findings, we conducted a detailed data 
analysis for ~15.25 yr Fermi-LAT observations on the y-ray 
excess and discovered a high-energy y-ray emission counter- 
part of a blazar, NVSS J171822+423948, exhibiting increas- 
ingly bright features in GeV 74-rays. Jiang et al. (2024) claimed 
that it is co-spatial with the IceCube neutrino event IC- 
201221A. This paper is structured as follows. In Section 2, we 
describe the Fermi-LAT data analysis procedures used in this 
work and report the y-ray results. We summarize and discuss 
some results in Section 3. 


2. Fermi-LAT Data Analysis 


In Zhang et al. (2023), a new y-ray source at 
R. A. = 17°18"29879 and decl. = +42°38/37”38 was found 
with a TS value of ~40, then it was cataloged as 4FGL J1718.5 
+4237 in 4FGL-DR4 based on Fermi-LAT 14 yr data (Ballet 
et al. 2023). We then carried out a data analysis for 4FGL 
J1718.5+4237 with Fermi-LAT 15.25 yr data. Using the tool 
gtfindsrc, we derived the position of 4FGL J1718.5+4237 at 
R. A. = 17°18™26388 and decl. = +42°38/30/46 with a 68% 
confidence level error radius of 2/1. The following analyses are 
based on this y-ray position. To identify the origin of the y-ray 
source, we carefully cross-checked all possible y-ray sources 
within 3/5 (2o0 error circle) around 4FGL J1718.5+4237 in 
Simbad database,? and totally retrieved eight objects including 
a blazar named NVSS J171822+423948 (Gaia Collaboration 
2020) with optical position of R.A.=17°18"22877 and 
decl. = +42°39'45”05 with angular distance of 1/5 from the 
y-ray position. We further investigated these sources with the 
NASA/IPAC Extragalactic Database (NED),* and excluded 
other known possible y-ray sources that are not found to be 
potential y-ray emitters. The nearest known y-ray source in the 
4FGL-DR4, 4FGL J1716.8+4310 (R.A. = 259°203, 
decl. = +43°173) associated with the y-ray counterpart of 
globular cluster NGC 6341, is about 0°58 away. Therefore, the 


3 https: //simbad.u-strasbg.fr/simbad/ 
4 http: //ned.ipac.caltech.edu / 
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most likely counterpart of the y-ray source 4FGL J1718.5 
+4237 would be NVSS J171822+423948, although we still 
cannot conclusively rule out the other sources within the 2c 
error circle of y-ray position. The following data analysis is 
aiming to determine the association between this y-ray source 
and NVSS J1718224+423948. 


2.1. Data Reduction and Best-fit Model 


The data utilized in the current study comprise all-sky- 
survey observations taken during 15.25 yr of Fermi-LAT 
operation, from 2008 August 4 (MJD 54682.6) to 2023 
November 1 (MJD 60249.3). The LAT Pass 8 Front + Back 
(evclass = 128 and evtype = 3) “Source” class events in the 
energy range 0.1-500.0 GeV were extracted from a 20° x 20° 
region of interest (RoI) centered on the y-ray position by the 
used Fermi Science Tools software package of Fermi- 
tools2.0.19. The events with zenith angles >90° were filtered 
out to avoid the y-ray contamination from the Earth’s limb. The 
expression “DATA_QUAL>0 && LAT_CONFIG = 1” was 
selected to obtain the good time intervals by employing the tool 
gtmktime. The instrument response function of “P8R3_SOUR- 
CE_V3” was adopted. 

We used a script named make4FGLxml.py° to generate a 
model file, which consists of y-ray sources within the Rol 
of the 4FGL-DR4 catalog. The model file included the 
spectral parameters, the spectral forms, and the coordinates 
for y-ray sources in the 4FGL-DR4 within 25° around 
4FGL J1718.5+4237. These initial values were derived from 
the catalog file of gll_psc_v34.fit.® The latest diffuse emission 
templates of gll_iem_v07.fits and iso_P8R3_SOUR- 
CE_V3_v1.txt were used to model the two y-ray backgrounds 
of the Galactic and extragalactic isotropic emissions,’ 
respectively. In the data analysis, the parameters of the flux 
normalizations and the spectral shapes were set free for the 
sources within 5°. While for the sources within 5°-10° or 
those outside 10° but identified as variable sources, we also 
freed their flux normalizations. The normalizations of the two 
backgrounds were set free too. We froze all other parameters 
at their values that were provided in 4FGL-DR4. For 4FGL 
J1718.5+4237, its y-ray spectrum is described by a power- 
law (PL) spectral shape of dN/dE = N(E/Ep)! with the 
spectral index [ = 2.5. 

Then we performed a binned maximum likelihood analysis 
between the whole database and the model file. The data 
reduction followed the standard criteria for the point-source 
analysis that is recommended by Fermi-LAT Collaboration.® 
The fit for 4FGL J1718.5+4237 results in TS ~ 135 in the 
0.1-500.0 GeV energy range, with an integrated average 
? https: //fermi.gsfc.nasa.gov /ssc/data/analysis /user/ 
© https: //fermi.gsfc.nasa.gov /ssc/data/access /lat/14yr_catalog/ 

7 https: //fermi.gsfc.nasa.gov /ssc/data/access /lat/BackgroundModels.html 
i https: //fermi.gsfc.nasa.gov /ssc/data/analysis /scitools/ 
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Figure 1. TS map (left) and residual TS map (right) with a 5° x 5° region around 4FGL J1718.5+4237 in 0.1-500.0 GeV obtained with ~15.25 yr of Fermi-LAT 
data. The coordinate of NVSS J171822+423948 is marked with a green cross. The y-ray source 4FGL J1718.5+4237 is tagged with a lime plus, and the lime circle 
stands for 2ø error circle of target. The position of 4FGL J1718.5+4237 reported in Jiang et al. (2024) is tagged with a purple plus (marked with “JX”). Other 4FGL- 
DR4 4-ray sources that begin with “4FGL” are shown with white diamonds. The TS value at each pixel on the sky is scaled with color. The two TS maps share the 
same color-bar for convenient comparison and have a spatial pixel size of 0°1 x O°1 on one side. 


Table 1 
Best-fit Results of the Likelihood Analysis for 4FGL J1718.5+4237 


Time Range (MJD) Model Parameter Value 

PL F TS Flux 
54682.6—60249.3 (Whole) 2.60(0.10) 135.09 6.80(1.30) 

2.50(0.12)* e e 

54682.6—57602.7 (Quiet) te 0.0 2.12° 
57602.7—60249.3 (Loud) 2.55(0.08) 226.25 11.92(1.72) 

LP a B E, (GeV) TS 

(Loud) 2.55(0.14) 0.13(0.07) 1.17 219.45 


Notes. The integrated photon flux is in the unit 107° photons cm~? 


à Parameter values of PL model provided in 4FGL-DR4, 
P The 95% (20) flux upper limit. 


photon flux of (6.80 + 1.30) x 107° photons cm~°s™' and a 
photon index of [~2.6. Its parameter values are in good 
agreement with those reported in 4FGL-DR4, except for the TS 
value. The best-fit results of the likelihood analysis for the 
target are listed in Table 1. The SEDs of AGNs in MeV-GeV 
ranges can also be described by a log-parabola (LP) model with 
a formula of dN/dE = No(E/E,) ¢+8s4/4») (e.g., Cerruti 
et al. 2013; Dermer et al. 2014, 2015). We also fitted the events 
with a spectrum of an LP model, and the best-fit results for the 
LP model are also shown in Table 1. Comparing the best-fit 
results of the two models, we believe that PL is the best model 


=j . 
s~ , and the errors are in parentheses. 


for the target. Then the best-fit results of the PL model were 
saved as a new model file. 

Based on the new model, we employed the tool gttsmap to 
calculate a 5° x 5° TS map centered at 4FGL J1718.5+4237 
in the energy range of 0.1-500.0GeV. In this step, we 
removed the target in the model file. All the parameters of the 
‘y-ray sources in the Rol, except the two diffuse components, 
were fixed at the values from the best-fit model. The TS map 
is displayed in the left panel of Figure 1. We mark the 4FGL 
J1718.5+4237 and NVSS J171822+423948 with a lime plus 
and a green cross, respectively. The lime circle stands for 20 
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Figure 2. The ~15.25 yr integrated flux of 4FGL J1718.5+4237 in 
0.1-500.0 GeV using a PL model with | yr time bins. The data points 
represent the source photon flux for time bins with TS > 5, and their TS values 
are displayed with pink bars. For TS < 5, we derived 95% flux upper limits, 
which are shown with gray arrows. The mean flux and the corresponding 
uncertainty range from the whole likelihood analysis are marked with a 
horizontal blue dash-dotted line and a blue shade, respectively. The fluxes and 
uncertainties of the Loud state are marked with a green dashed line and a green 
shade, respectively. Note that the last data point only covers 0.25 yr. 


error circle of the y-ray position of 4FGL J1718.5+4237. In 
order to examine contamination from the nearby potential 
y-ray source (not included in 4FGL-DR4), a residual TS 
map was obtained based on the above best-fit model. The 
residual TS map is featured in the right panel of Figure 1. We 
found that no pixel having TS > 25 is seen in the residual 
TS map. 


2.2. Long-term Flux Variability 


To investigate the temporal variability behavior in GeV, we 
constructed a light curve for the target in 0.1-500.0 GeV by 
performing an unbinned likelihood analysis based on the new 
model file. In order to show the light curve as completely as 
possible and to facilitate future comparisons between our data 
and others’ results, the yearly binning method is used in this 
step. Here, we fixed the parameters of spectral shapes at their 
best-fit values and only freed the normalizations mentioned 
above for the sources in the RoI. We presented the flux and the 
corresponding TS value for each time bin in Figure 2. The data 
points with TS > 5 were included to show the light curve as 
completely as possible. Otherwise we derived their 95% flux 
upper limits for the data points with TS < 5 and colored them 
gray. Additionally, it should be noted that the last data point is 
incomplete, covering only 3 months. From Figure 2, we know 
that the target’s y-ray emissions are mainly concentrated in the 
latter part, reaching its maximum -ray flux state in the last 
time bin recorded so far. 
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2.3. Quiet and Loud States 


According to Figure 2, we defined the time period as two 
distinctive emission states: Quiet and Loud states, marked with 
gray and black shades in Figure 2, respectively. A Loud state 
refers to the period of detection of high flux from 2016 August 
(MJD 57602.69) to the end of the data, and a Quiet state is 
selected as the time period when the source is in a relatively 
fainter state from 2008 August to 2016 August. We also carried 
out the likelihood analysis for the two states. For the Quiet state, 
we obtained the TS value of 0 and a flux upper limit of 
~2.12 x 107° photonscm™*s~', and the Loud state results in 
TS = 226, and the flux is ~1.19 x 1078 photons cm~? s~! with 
T ~ 2.6 in 0.1-500.0 GeV. The results are presented in Table 1. 
The index T is consistent with them in the whole fit result, which 
suggests that no spectral variation occurs during the Loud state. 
All the best-fit results were saved as two updated model files, 
which will be used in the next procedure of the analysis. 

Based on two updated models, we created 5° x 5° TS maps 
for the two states, respectively. They are shown in the left and 
right panels of Figure 3. These two TS maps used the one color 
bar to scale their TS values for each grid location on the sky. 


2.4. Spectral Analysis in Loud State 


To investigate the spectral shape in 0.1—500.0 GeV during the 
Loud state, we constructed its SED based on the updated model. 
We fixed the parameters of spectral shape for all the sources at 
their best-fit values and only freed the normalizations of the 
target and two backgrounds. We divided the 0.1-500.0 GeV 
energy band into 12 equally logarithmically spaced energy bins, 
and then performed the likelihood analysis to calculate the flux 
and TS values for each energy bin. The y-ray data points with 
TS > 5 are represented in Figure 4 with black marks, otherwise 
95% upper limits of flux are shown with gray marks. The PL 
best-fit model is displayed with a blue solid line. 

The LP model was also used to fit the y-ray events from the 
target. We obtain a spectral slope a ~ 2.6, a curvature parameter 
around the peak 3~ 0.13, and a TS ~ 220 (listed in Table 1). 
We show them in Figure 4 with an orange dash-dotted line. A 
likelihood ratio test was used to check the PL model (null 
hypothesis) against the LP model (alternative hypothesis). We 
compared these values by defining the curvature test statistic 
TScurve = 2log(Lip/Lp_) (Nolan et al. 2012), and the TSeurve 
was calculated to be ~ — 3. This maybe indicates no statistical 
evidence for a curved spectral shape. From Figure 4, we also can 
see that the PL model gives an acceptable description for the 
SED, which agrees well with the whole fitting. 


3. Results and Discussion 


In this work, we present the preliminary results of our study 
on 4FGL J1718.5+4237 in the latest release of LAT source 
catalog 4FGL-DR4 (Ballet et al. 2023) using the Fermi-LAT 
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Figure 3. Two TS maps in the left and right panels are obtained with the events from Quiet and Loud states, respectively. Other aspects are the same as in Figure 1. 
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Figure 4. 0.1-500.0 GeV SED of 4FGL J1718.5+4237 in the Loud state. The 
95% flux upper limit is calculated for the data points with TS < 5. The blue 
solid and orange dash-dotted lines represent the best-fit results of PL and LP 
models, respectively. 


Pass 8 data set covering 15.25 yr from 2008 August to 2023 
November in the energy range 0.1-500.0 GeV, and it was 
detected at a significance level of ~11.6c. This y-ray source is 
positionally associated with the blazar NVSS J171822 
+423948 at a high-redshift of ~2.7. Our results suggested 
that NVSS J171822+423948 is a new y-ray emitter and has 
remained in quiescence for a long time. In the derived TS maps 
(see Figure 1), the position of NVSS J171822+423948 falls 


into the lo error circle of 4FGL J1718.5+4237 in y-rays. The 
angular separation between the optical and y-ray coordinates is 
derived as 1/5. This observation provides further evidence of 
NVSS J171822+423948 being a new y-ray emitter in the 
extragalactic high-energy sky. 

From the long-term light curve, interestingly, the source 
showed a significant increase in y-ray flux and continuous 
brightening since approximately MJD 57602.69 (date 2016 
August 2 16:29:16.8 UTC). We performed a spectral analysis 
of the y-ray emissions from NVSS J171822+423948 during 
the Loud state. The results are listed in Table 1 and the SED is 
displayed in Figure 4. The y-ray spectrum is well described by 
a PL model. The y-ray flux in the brightest period (the last data 
point) increased as much as ~12 orders of magnitude 
compared to the flux upper limit in the Quiet state (estimated 
over 8 yr of Fermi-LAT observations). Its y-ray spectrum is 
soft, which aligned with that typically observed from high- 
redshift blazars (z > 2) (e.g., Ackermann et al. 2015). 

In the previous literatures (Véron-Cetty & Véron 2006, 
2010; Ahn et al. 2012; Gaia Collaboration 2020), they 
presented a compilation of all known AGNs, and NVSS 
J171822+423948 was one entry in the catalog and classified as 
a blazar. In Zhang et al. (2023), a significant y-ray excess with 
a TS of ~40 was presented in their results, but they did not 
provide much discussion on this source. It was subsequently 
cataloged in 4FGL-DR4 (Ballet et al. 2023). Jiang et al. (2024), 
reported on the multi-wavelength and multi-messenger studies 
of this blazar, including its y-ray emissions, and claimed that 
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this blazar is the low-energy counterpart of 4FGL J1718.5 
+4237. They also stated that 4FGL J1718.5+4237 is the only 
gamma-ray source known to be co-spatial with the IceCube 
neutrino event IC-201221A with an arrival time at MJD 59204, 
which corresponds to 2020 December 21. This will be very 
important for the study of high-energy neutrino astronomy. 
Their analysis of LAT data extends up to MJD 59945 (2023 
January 1). Moreover, this y-ray source exhibited no variation 
in its GeV light curves from MJD 58856 (2020 January 8) to 
59945 (2023 January 1). Here, we extend the time range of the 
data analysis and found a trend of gradually increasing gamma- 
ray emissions from 4FGL J1718.5+4237 starting from MJD 
58333 (2018 August 2). Therefore, further multi-wavelength 
observing campaigns for NVSS J171822+423948 are strongly 
encouraged, which will be important to reveal the properties of 
this high-redshift blazar. Particularly, carrying out optical and 
X-ray observations of this y-ray source during the Loud state 
will enable us to compare the broadband SED and constrain the 
underlying astrophysical mechanisms. 
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